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In a spirit akin to the sandpile model of self-organized criticality, we present a simple statis- 
tical model of the cellular-automaton type which produces an avalanche spectrum similar to the 
characteristic-earthquake behavior of some seismic faults. This model, that has no parameter, is 
amenable to an algebraic description as a Markov Chain. This possibility illuminates some impor- 
tant results, obtained by Monte Carlo simulations, such as the earthquake size-frequency relation 
and the recurrence time of the characteristic earthquake. 

I. INTRODUCTION 

If there is a well-established fact about regional seismicity this is the relationship between the magnitude of an 
earthquake and its frequency, known as the Gutenberg- Richter (GR) law jy . This law is of the power-law type when 
magnitudes are expressed in terms of rupture area 

NotS- b , (1.1) 

where N is the number of observed earthquakes with rupture area greater than S, and b is the so-called 6-value, 
which is a "universal constant" in the range 0.8-1.2 J2|. The GR law implies that earthquakes are, on a regional or 
world-wide scale, a self-similar phenomenon lacking a characteristic scale (but see ||). 

It is important to notice, however, that the GR law is a property of regional seismicity, appearing when we av- 
erage seismicity over big enough areas and long enough time intervals. Recently, a wealth of new data has been 
collected to extract statistics on individual systems of earthquake faults Q. Interestingly, it has been found that 
the distribution of earthquake magnitudes may vary substantially from one fault to another and that, in general, 
this type of size-frequency (SF) relationship is different from the GR law. Many single faults or fault zones display 
power-law distributions only for small events, which occur in the intervals between roughly quasi-periodic earthquakes 
of much larger size which rupture the entire fault. These earthquakes are termed "characteristic", and the resulting 
SF relationship, characteristic earthquake (CE) distribution. 

There is much debate about the origin of the CE distribution ||. Because of the short period of instrumental 
earthquake records and the scarcity of paleoseismic studies Q], the statistics of naturally occurring earthquakes in 
single faults are poor. This fact justifies the development of "synthetic seismicity" models H, in which long catalogs 
of events are generated by computer models of seismogenesis. Such models can be tuned by requiring that they 
reproduce what is known of the statistics of past seismicity to a reasonable degree, and then use them to forecast 
statistical inferences about the behavior of seismicity using much longer and homogeneous catalogs of synthetic events. 

Many different seismicity models have been presented in the past twenty years or so. Robinson and Benites Q| 
classify these modeling approaches into five groups: (1) cellular automata models, (2) spring-block models, (3) models 
of single faults in which slip is discretized into patches and obey simplified friction laws, (4) continuum models that 
I* ■ utilize realistic constitutive friction laws, and (5) actual physical models. 
. £h ! Cellular automata models H have only recently appeared in seismological literature, hand in hand with the concept 
of self-organized criticality These models are usually nondeterministic and represent faults as one- or two- 
dimensional features. They neglect the details of both elasticity and fault friction, substituting them by simple 
cellular automata rules. Despite their simplicity, they are able to reproduce various types of SF statistics, including 
GR and CE distributions @. 

The key ingredients of any of these models are: (1) the dimensionality of the fault (ID or 2D), (2) the number of faults 
included in the model (one, a few, or many faults), (3) the employed stress transfer mechanism (nearest-neighbors, 
long-range elasticity, mean-field), (4) the degree of incorporation of inertial effects (quasi-static, quasi-dynamic, or 
fully dynamic), (5) the assumed constitutive stress-slip law (experimental, static-dynamic, velocity-weakening, etc.), 
and (6) the degree of stress conservation (conservative versus dissipative models). 

Various discrete models of seismicity are able to show a transition from GR to CE statistics. Among them, we 
want to cite here the models of Ceva and Perazzo [Ol, Carlson et al. I|13l, Lomnitz-Adler |14|, Rundlc and Klein 
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Jig ], Ben-Zion and Rice [jL5||, Dahmen et al. ||, Moreno e£ a/. |T^] and Heinz and Zoller jig]. All of them have in 
common a transition from a GR to a CE distribution when some of the model parameters are continuously changed. 



17LP.81 ); others appeal to percolation theory ( |_2[); 
16[). But there is no consensus about the origin of 



Some authors interpret the results in the light of SOC (e.g. []13l 
and still others talk about generic critical point behaviour ( [nj 
these spatiotemporal complexities, with two extreme views: the view that says that the spatiotemporal complexity is 
generated by the dynamics of the fault, and the opposite view that favors nonuniformitics in geometry and material 
properties as the cause of this complexity. 

Our purpose here is to build the simplest cellular automaton model of seismicity capable of displaying a SF rela- 
tionship of the CE type. That is, a model which exhibits a power-law relationship for small events and an excess 
of big event (of the order of the system size), together with a very low probability of events of intermediate size. 
With respect to the six basic ingredients of discrete models of seismicity introduced above, the model presented here 
is (1) one-dimensional, (2) for a single fault, (3) with a percolation-like stress-transfer mechanism, (4) quasi-static, 
(5) static/dynamic with total stress drop, and (6) dissipative. Also, because of the inherent simplicity of the model, 
we want to be able to derive analytically some of the statistical properties of the resulting synthetic seismicity using 
Markov chains. 



II. THE MODEL AND ITS SIMULATIONS 



Consider a one dimensional vertical array of length N . The ordered positions, or levels, in the array will be labeled 
by an integer index i varying from 1 to N. This system performs two functions, it is loaded by receiving individual 
particles in the various positions of the array, and unloaded by emitting groups of particles through the first level, 
i = 1, which are called avalanches or earthquakes. 

These two functions proceed with the following four rules: 

(i) The incoming particles arrive at the system at a constant time rate. Thus, the time interval between each two 
succesive particles will be considered the basic time unit in the evolution of the system. 

(ii) All the positions in the array, from i = 1 to i = N, have the same probability of receiving a new particle. When 
a position receives a particle we say that it is occupied. 

(iii) If a new particle comes to a level which is already occupied, this particle disappears from the system, or in other 
words, this particle assignment is wasted. Thus, a given position i can only be either non-occupied when no 
particle has come to it, or occupied when one or more particles have come to it. 

(iv) The level i = 1 is special. When a particle goes to this first position an avalanche occurs. Then, if all the 
successive levels from i = 1 up to i = k are occupied, and the position k + 1 is unoccupied, the effect of the 
avalanche is to unload all the levels from i = 1 up to i = k. Hence, the size of this avalanche is k, and the 
remaining levels i > k remain unaltered in their occupancy. 

Thus, from what has been mentioned above, this model has no parameter; the size N is the unique specification to 
be made, and the spatial correlation is induced by the ivth rule. Now, the state of the system is given by stating which 
of the (i > 1) N — 1 levels are occupied. Each of these states corresponds to a stable configuration, and therefore the 
total number of possible configurations is 2^ N ~ X \ We use the term "total occupancy" for the configuration in which 
all but the first level are occupied. 

After the occurrence of an avalanche the system is left in a stable configuration; the following particle additions go 
progressively loading the system, and when a particle is again assigned to the first level a new avalanche is triggered. 
Each avalanche empties the lower levels of the system as explained in rule (iv), and the system is left in another stable 
configuration. The size of the avalanches can range from 1 up to N and the avalanche of maximum size, k — N, will 
be called the characteristic one. 

From these evolution rules we deduce that after a time unit, i.e., after a new incoming particle assignment, we 
will have an avalanche if the new particle goes to i = 1, and this occurs with a probability 1/N. Conversely, with 
a probability of (N — 1)/N there will be no avalanche. In this case the system will advance one unit in its level of 
occupation when the new particle is assigned to a non-occupied level, and it will remain at the same configuration if 
the assigned level was already occupied. 

As this model is 1-dimensional, extensive Monte Carlo simulations can be performed to accurately explore its 
properties. In this paper we focus on two important properties, the avalanche size-frequency relation and the statistics 
of the time of return of the maximum-size avalanche. 
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The results for the avalanche size-frequency relation, pk, are drawn in Fig. [l] and written in Table M. In Fig. |lj we 
have superposed pk for TV = 10, TV — 100 and TV — 1000 . This has been partly included in Table | as well. In Fig. |l| 
there are three notable properties to be commented on. First and most important, we see that the characteristic 
avalanche, k = TV , has a much higher probability of occurrence than the avalanches of big size but with k < TV. In 
fact, for TV = 10, 100 and 1000, the probability of their respective characteristic avalanches does not differ much, and is 
near 10%. We can express this fact by saying that in this model, grosso modo, one would likely only observe very small 
avalanches and the characteristic one. Secondly, forgetting for a moment the case k — N , we observe an approximate 
power law behavior, a la Gutenberg-Richter, for the rest of avalanches. The exponent b of this differential distribution 
is roughly 1.6. And thirdly, we observe the perfect coincidence of these curves of probability for systems of different 
size TV. This is also appreciated in the numbers collected in Table |[ This, in a sense, unexpected size-invariance will 
be discussed in detail in the next Section. 

In Fig. U we represent the probability curve for "the time of recurrence of the characteristic avalanche" . This curve, 
obtained by Monte Carlo simulations, corresponds to a system of TV = 10. In the abscissas axis, time (denoted by n) 
starts at just after the occurrence of a k — TV avalanche. It is clear in this model that only after a minimum time 
lapse of size TV the probability of occurrence of a new k = TV avalanche can be non null. We observe in Fig. |^ that 
after this minimum time lapse, P(n) grows to a maximum and then decays. For the size TV = 10 analyzed in this 
figure, the maximum of probability corresponds to a time interval n = 34. 



III. THE MODEL AS A MARKOV CHAIN 



It is easy to become convinced that, for a given TV, the 2( Ar_1 ) stable configurations of our model can be considered 
as the states of a finite, irreducible and aperiodic Markov chain with a unique stationary distribution jigj] . These 
configurations are classified in groups according to its occupation number (number of occupied levels); the number 

of configurations with j occupied levels is C ( ^ ^ ^ . One step in the chain corresponds to the result of adding 

a new particle to the system. Up to approximately TV = 10, the transition matrix, M, can be easily obtained using 
Mathematica as well as the corresponding stationary probabilities for each configuration, which correspond to the 
components of the eigenvector of M with eigenvalue, A, equal to unity. 

For small TV, M and its eigenvectors are obtained by inspection. Let us then start reproducing the first numbers 
quoted in Table || for the probabilities of occurrence of avalanches in systems of small size. With this aim, let us 
consider Fig. a. There, in (A), (B) and (C) there appear all the stable configurations, ordered in an increasing state 
of occupation, for TV = 2, TV — 3, and TV — 4, respectively. For the moment, the black level in the top position of the 
configurations has no meaning. 

For TV = 2, using the same order and notation for the configurations as in Fig. ||A, the transition probabilities are 
M 0j0 = 1/2, M a , b = 1/2 , M M = 1/2, and M bfi = 1/2. Thus 



M 



1/2 1/2 
1/2 1/2 



(TV 



This M is a doubly stochastic matrix and hence the two stationary probabilities are equal. 

p. = 1/2, Pi = 1/2, (TV = 2). 



(3.1) 



(3.2) 



For TV = 3, the non-null transition probabilities are: M a ^ a — 1/3, M a ^ — 1/3, M a c = 1/3; Af&.b = 2/3, M},,d = 1/3; 
M c , = 1/3, M c , c = 1/3, M c4 = 1/3; M d . a = 1/3, and M dA = 2/3. Thus 



M 



1/3 1/3 1/3 \ 

2/3 1/3 

1/3 1/3 1/3 

1/3 2/3 / 



(TV = 3). 



And the components of the eigenvector corresponding to A = 1 are: 

Pa = 1/4, p 6 = l/8, p c = l/4, p d = 3/8, 



(TV = 3). 



(3.3) 



(3.4) 



Finally, for N = 4 the non- null transition probabilities are M a a = 1/4, M a 6 = 1/4, M a c — 1/4, M a d — 1/4; Mj & 
2/4,Af b , e = 1/4, M bJ = 1/4;M C , C = 2/4, M c , e = l/4,M c>fl ='l/4,M d , = l/4,M d , d = l/4,M dJ = 1/4, Ai^ 
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1/4; Af e , e = 3/4, M eJl = 1/4; M f , b = 1/4, M fJ = 2/4, M /)h = 1/4; M g , a = 1/4, M 5iS = 2/4, M g , h = 1/4; M h . a = 
1/4, M M = 3/4. Thus 



M = 
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(N = 4). 



(3.5) 



And, after its diagonalization, one finds 



Pa = 9/64, p fo = 7/64, p c = 9/128, Pd = 3/64, 
p e = 23/128, p/ = 5/64, p 3 = 15/256, p h = 81/256, (iV = 4). 



(3.6) 



From these numbers, one deduces that in a system with N = 2 levels one should expect avalanches of size k = 1 with 
a probability pi = p a = 1/2 , and of size k = 2 with P2 = Pb = 1/2. In iV = 3, pi = p a +Pb = 1/2, P2 = Pc = 1/8, and 
p 3 = Pd = 3/8. And in N = 4, pi = p a + p b + p c + p e = 1/2, p 2 = Pd+Pf = 1/8, p 3 = p g = 15/256 = 0.05859.., and 
Pi = p h = 81/256. 

Thus, we have observed that in systems with N — 2, N = 3 and N = 4 levels, the value of pi is a constant equal 
to 1/2 and this result coincides with Table |[ We have also observed that in N = 3 and N = 4 the value of P2 is a 
constant equal to 1/8, and in N = 4 we have deduced that P3 = 0.05859. All this agrees with Table |. 

A conclusive argument proving that in this model the value of pk is a constant independent on the size N is obtained 
by re-analyzing Fig. ^| from a wider perspective. Now we consider that in Fig. |j|A, the configuration labeled by a 
represents, in a system of size N , all the configurations in which the level i = 2 is not-occupied; the rest of the levels 
i > 2 are in any possible state of occupation (this is represented by the black level at the top of the configuration). 
And in Fig. ||A, configuration b, we consider a system of size N which has its second level occupied. The probabilities 
of these two cases must add to unity. Then defining a Markov chain for these two excluding states, and using the same 
notation as before, we find M a ^ a = (N - 1)/N, M a , 6 = 1/iV, M _ a = 1/N , and M M = (N - l)/N. The diagonalization 
of this matrix leads for p a and p b to the same results written in Eq. (3^3), p a = Pb = 1/2; however, its interpretation 
now is different. Here, p a — 1/2 implies that for any value of N, the probability of having an avalanche k = 1 is 1/2; 
and pt, = 1/2 simply means that the probability of having avalanches k > 1 is 1/2. Using the same line of reasoning, 
and referring to Fig. configuration a represents all the configurations, in a system of size N, where levels i = 2 
and « = 3 are not occupied. And configuration b represents all the configurations where the level i = 2 is free and the 
level i = 3 is occupied, etc. Then, the non null transition probabilities are M a a = (N — 2)/N, M a {, = 1/N, M a c = 
1/N; M b , b = (N- 1)/N, M M = 1/N; M c , a = 1/N, Af c , c = (N - 2)/N, M cA = 1/N; M d , a = 1/N, M dA = (N - 1)/N. 
The diagonalization of this matrix provides the same stationary probabilities quoted in Eq. ( |3.3[ ). Here p c = 1/8 
means that for an arbitrary N, the probability of occurrence of avalanches of size 2, pi is 1/8. And the fact that 
Pa + Pb = 1/2 confirms that for any N, p\ = 1/2. 

Extending this line of reasoning to the 8 configurations drawn in Fig. ^C, one concludes that, for any N, P3 = pg — 
15/256 = .058593.. and one verifies the previous conclusions pi = 1/2 and P2 = 1/8. 

Therefore, in this model, if for the system of size N one knows all the pk from k = 1 to k = N, then for the system 
of size N + 1 the pk are identical , with the exception of the last two, and these fulfill 



p N {N) = PN+1 (N + 1) +p N (N + 1) 



(3.7) 



the recursive way in which ppf(N) divides into pn+i(N + 1) and pn(N + 1) is, however, non trivial. 

Let us analyze now, from the Markov-chains point of view, the results for the time of return of the characteristic 
earthquake shown in Fig. |^. After a k = N avalanche, the system is left in the configuration of no occupancy (for the 
present discussion we will refer to this configuration as a{). A new characteristic avalanche will occur when, starting 
from the configuration ai, the system reaches the configuration of total occupancy (which will henceforth be denoted 
by a at) j and then the next particle is assigned to the i = 1 level. The number of time steps elapsed between a\ and 
apf, plus 1, will be denoted by n. And our purpose is to compute the probability of occurrence of a k = N avalanche 
as a function of n. It must be understood that between a\ and the occurrence of the next k = N avalanche on time 
n, the system may have visited ojv an arbitrary number of times but without triggering any k — N avalanche. In 
other words, in those visits there has been no transition from to a\. 
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In Markov chains, the transition matrix M gives the probability of going from one configuration to another in one 
step, and the m step transition probability is the m-th power of M. Thus a simple way to compute P(n) is the 
following: 

1 Take M, point to the element in the last row and the first column, and substitute it by 0. We will denote the 
new matrix M' . 

2 Compute T n = M l[n ~ 1] . 

3 Take the element of the first row and the last column of T n , and multiply it by 1/N. This result is P(n). 

This is clear, because M' does not permit transitions from ajv to oi. Thus, in T n we have all the probabilities 
of transitions, in n — 1 steps, between all the configurations, with the restriction that from to ai this transit is 
forbidden. Hence, in T n the matrix element of the first row and the last column corresponds to the transition from 
a% to ajy, in n — 1 steps and with the mode ajv — ► ai locked. Finally, as the 1/N factor is the probability of the 
transition ajv — ► ffli, we actually have built P(n). In Fig. ||, for N — 10, we see the perfect match between the Monte 
Carlo simulations and the results coming the theory of Markov, calculated using Mathematica. 

In ord er to get a quantitative insight on P(n), let us apply this method to N = 2. In this case, M is given by 
Eq. (|3J1) . Hence, 

and, therefore, 

P{n) = (n - 1) • (l/2)("-« ■ (1/2) = (N = 2). (3.9) 

Thus, we observe that the asymptotic fall-off behavior of P(n) is of the type 

P(t) octexp(-i), (N = 2). (3.10) 

It is important to recall that in aperiodic, irreducible and finite Markov Chains such as that of our model, the mean 
waiting time for a configuration is the inverse of the stationary probability of that configuration. Then, the mean 
time between characteristic avalanches in this model is 

IN , 

<">= — = —■ (3- 11 ) 
— ojv 

As an example, for N — 4, < n >= 4 ■ (256/81) = 12.4 time units. 



IV. CONCLUSIONS 



We have presented a one-dimensional discrete model of seismicity that displays a size-frequency spectrum similar to 
that expected from a characteristic-earthquake behavior. Although this stochastic model obviously lacks the explicit 
physics of other detailed -and complex- dynamical models of earthquakes, its basic hypotheses and implications are 
clear, phenomenologically reasonable and coherent. The size invariance of this model, also surprisingly manifests itself 
in predicting an identical probability of occurrence for earthquakes of the same magnitude independently of the size 
N of the system. In this universal rule the characteristic earthquake is excluded. This model has the additional bonus 
that several important predictions can be algebraically derived by using the theoretical framework of the Markov 
chains. Specifically, the statistics of the time of return of the characteristic earthquake are neatly predicted by this 
formalism. 

Dahmen et al. || report a model able to transit from the GR to the CE behavior and, what is more interesting, 
these authors define a "configurational entropy" as an appropriate concept that reflects in which of these two extreme 
behaviors the system is actually operating. In qualitative terms, the GR behavior corresponds to a high entropy 
mode of operation while the CE behavior corresponds to a low entropy mode. Therefore, we would like to make an 
assesment of our model from this configurational entropy point of view and check the concordance, or not, with the 
conclusions of Ref . || . 
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In our model, the configurations are classified in groups according to the number of levels, j, that are occupied. 

The statistical weight of each j is C f ^ J, which has its maximum values for j around N/2. Conversely, the 

statistical weight is minimum on the extrema: for j around 1 and for j near N — 1. Then, we need to find out 
the values of the stationary probabilities of each occupation number j. This is shown in Fig. for a system of size 
N = 100. There we observe how in our model the system resides most of the time in the configurations of maximum 
occupancy, that is, where the configurational entropy is a minimum, agreeing with the interpretation of Ref. [[|. 

As a final minor remark, note that in Fig. ^ P(j = n) and P(j = N-i) are identical. This is a property that holds in 
this model for any value of N and which can be easily proved. For brevity reasons, we omit this proof. 
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FIG. 1. Probabiliy of occurrence of earthquakes of magnitude k. The results for N — 10, 100, 1000 are superposed. 
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FIG. 2. Probability of return of the characteristic earthquake as a function of time. 
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FIG. 3. Explicit configurations for: A)N = 2, B)N = 3 and C)N = 4. 
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FIG. 4. Stationary probabilities as a function of the number of levels occupied, j. 
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TABLE I. Probability of occurrence of the earthquake of magnitude k for a different system size, N. 
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